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We  consider  the  problem  of  reconstructing  an  image  from  a  noisy  record.  We 
describe  existing  methods  due  to  Geman  and  Geman  (1984)  and  Besag  (1986)  which  use 
a  Markov  random  field  model  for  the  true  scene  but  assume  that  each  pixel  consists  of  a 
single  colour.  In  order  to  improve  the  quality  of  the  restoration  at  the  boundary  of 
regions  of  different  colours  we  extend  these  methods  to  allow  pixels  to  contain  two 
regions  of  colour  separated  by  a  single  straight  line.  An  algorithm  for  performing  the 
reconstruction  is  presented  and  illustrated  by  an  example.  ^  j  ' 

l.  Introduction 


We  consider  a  rectangular  region  partitioned  into  pixels  labelled  1,2 . n.  Each 

pixel  is  coloured  black  or  white  and  the  colour  of  pixel  i  is  denoted  by  jq  which  takes 
the  value  0  for  white  and  1  for  black.  The  xi  are  unobserved.  We  work  instead  from  the 
observed  record  y,  which  consists  of  x,  plus  added  noise.  We  denote  the  whole  scene  by 
x  =  {jq;  i‘=l,...,n}  and  the  set  of  records  by  y  -  {y,-;  i=l,...,n}.  The  noise 

distribution  will  be  assumed  to  be  known  but  if  this  were  not  the  case,  it  could  be 

established  by  studying  training  data. 

Recent  developments  in  statistical  restoration  methods  use  a  Bayesian  approach. 
The  maximum  a  posteriori  (MAP)  estimate  of  the  true  scene  is  the  value  of  x  which 

maximises  Pfxly),  the  conditional  probability  of  x  given  the  record  y.  By  Bayes’ 

theorem 


P(x|y)  oc  /(y|x)  p(x ), 


(1.1) 


where  /(y|x)  is  the  conditional  likelihood  of  the  observed  record,  y,  given  the  true 
colouring,  x,  and  p(x)  is  the  prior  probability  of  x. 


DUmiBUTION  STATEMENT  A 

Approved  ter  public  release! 
Distribution  Unlimited 


We  assume  the  conditional  density  function  /(y(|x,)  to  be  known  and  for  the 
remainder  of  this  paper  we  shall  assume  that  the  records,  yt,  are  independently 
distributed  as  Gaussian  with  mean  xi  and  variance  a2.  Thus, 


Ky\x)  =  n/O’.W  =  (2™2)  2 exPi—\'L^i~xi)2)- 

1=1  2cr  i=  i 

-1  'A 

To  obtain  a  valid  formula  for  p(x),  we  assume  that  the  true  scene  corresponds  to  a 
locally  dependent  Markov  random  field  (MRF)  with  respect  to  a  specified  neighbourhood 
system,  that  is,  the  conditional  distribution  of  pixel  i  given  the  colourings  of  all  other 
pixels  depends  only  on  the  neighbours  of  pixel  i.  We  shall  use  a  second  order 
neighbourhood  system  in  which  pixels  are  considered  to  be  neighbours  if  they  are 
horizontally,  vertically  or  diagonally  adjacent  to  each  other.  A  detailed  definition  and 
further  examples  of  Markov  random  fields  may  be  found  in  Besag  (1974). 

The  form  of  p(x)  is  determined  by  the  nature  of  the  Markov  random  field.  In  our 
case,  we  have 


p(x)  «  e 

where  Z(x)  is  the  number  of  discrepant  pairs  in  the  scene,  x,  i.e.  the  number  of  pairs  of 
neighbours  which  are  of  opposite  colour,  and  P  is  a  fixed  positive  constant  (normally 
chosen  to  be  between  0.5  and  1.5). 
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The  maximisation  of  P(x\y)  now  corresponds  to  the  minimisation  of 

-ViO',-*,)2  +  PZOO  (1-2) 

2<r-i=i 

over  values  of  x  =  {xt;/=l 

This  expression  may  be  regarded  as  a  penalty,  the  first  term  penalising  any 
difference  between  the  record  and  the  fitted  value,  the  second  term  penalising  excessive 
roughness  in  the  reconstruction.  Clearly,  with  2n  possible  values  for  x  this  is  a 
computationally  large  problem  and  necessitates  the  use  of  a  sophisticated  algorithm. 

Geman  and  Geman  (1984)  use  the  method  of  simulated  annealing  which  attempts  to 
find  the  MAP  estimate  of  x  given  the  record  y.  Their  method  is  computationally 
extravagant  and  more  recent  developments  by  Greig,  Porteous  and  Seheult  (1986)  show 
that  the  MAP  estimate  of  any  two  colour  scene  may  be  found  exactly  using  the  Ford- 
Fulkerson  labelling  algorithm  for  maximising  flow  through  a  network. 
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Besag  (1986)  proposed  the  computationally  simpler  method  of  iterated  conditional 
modes  (ICM)  which  updates  each  pixel  in  turn,  choosing  for  it  the  most  likely  colour 
based  on  its  record  and  the  current  colouring  of  its  neighbours.  In  updating  pixel  i  the 
new  x ,  is  chosen  to  minimise  the  sum  of  terms  involving  xi  in  the  penalty  (1.2),  i.e. 
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■(yi-x,)2  +  pz(x,) 


where  Z(x()  is  the  number  of  neighbours  of  pixel  i  in  the  current  restoration  which  are 
of  the  opposite  colour  to  x,.  The  method  proceeds  by  scanning  the  scene,  successively 
updating  each  pixel  until  convergence  is  reached.  This  will  normally  occur  at  a  local 
rather  than  global  maximum  of  P(x|y),  but,  given  the  possibility  of  undesirable  long 
range  dependencies  in  the  MRF  model,  this  is  not  a  serious  drawback  and  might  even  be 
an  advantage. 


2.  Split  Pixels 

So  far  we  have  considered  scenes  in  which  each  pixel  is  coloured  wholly  one 
colour.  We  now  allow  pixels  in  the  true  scene  to  be  coloured  partly  black  and  partly 
white.  Each  record  y,  is  distributed  as  Gaussian  with  variance  o2  and  mean  pit  the 
proportion  of  pixel  i  which  is  coloured  black.  The  restoration  methods  that  we  have 
previously  discussed  can  be  used  for  this  problem  by  proceeding  as  if  the  pixels  were 
only  of  one  colour  but  the  quality  of  the  restoration  at  the  edges  of  objects  or  regions 
will  obviously  be  poor.  Instead,  we  can  allow  pixels  in  the  restored  image  to  be 
coloured  partly  black  and  partly  white.  The  simplest  form  of  this  is  to  quarter  each  pixel 
and  allow  it  to  be  filled  with  th,*.  most  likely  of  the  24  configurations.  This  method, 
proposed  by  Jennison  (1986)  uses  a  modified  version  of  ICM,  firstly  iterating  at  full 
pixel  size  and  subsequently  restoring  the  quarters;  in  the  second  stage  the  same  form  of 
MRF  model  is  used  for  the  subpixels  as  is  originally  used  for  full  pixels  This  method 
appears  to  work  well  and  has  prompted  work  into  the  further  breakdown  of  pixels. 

For  further  refinement  we  can  either  (i)  consider  an  mxm  breakdown  of  each  pixel 
or  (ii)  use  continuous  lines  within  the  pixel  to  represent  the  edge.  The  implementation 
of  (i)  requires  the  minimisation  of 

in  im#»  It  m  m 

— riO'.  +  z(xijk), 

l=  i  m  j—  i  k-  [  /=  i  j=  i  t=i 

where  the  subscript  ijk  refers  to  subpixel  j,k  within  pixel  i;  Xyk  takes  value  Owl  and 
Z(x,y* )  is  the  number  of  subpixel  neighbours  of  subpixel  ijk  in  the  current  restoration 
which  are  of  the  opposite  colour  to  x^k.  Note  that  subpixels  at  the  edge  of  a  pixel  will 
have  some  subpixel  neighbours  contained  in  an  adjacent  pixel.  We  can  see  that  as  m 
increases  this  minimisation  becomes  computationally  cumbersome.  Also,  it  offers  only  an 
approximation  to  (ii)  and  it  turns  out  to  be  easier  to  pass  to  the  limit  and  work  directly 
with  continuous  solutions. 

The  most  basic  form  of  (ii)  allows  a  single  straight  line  edge  within  each  pixel  and 
it  is  the  implementation  of  this  that  we  shall  describe.  It  is  no  longer  meaningful  to  talk 
of  discrepant  pixel  or  subpixel  pairs  and  we  replace  the  second  term  of  (1.2)  by  a 
multiple  of  the  total  length  of  edge  in  the  reconstruction  x.  Thus,  the  restored  image  is 
chosen  to  minimise 
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— T  £0’, -p,(*))  +  p  'L(X), 
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over  images  x  made  up  of  pixels  xt ,  i=l,...,n  either  of  a  single  colour  or  divided  into 
two  regions  of  different  colours  by  a  single  straight  line;  p,( jc)  denotes  the  proportion  of 
black  in  pixel  i;  L(x)  is  the  total  edge  length  in  scene  x  and  p'  is  a  fixed  constant  related 
to  the  p  used  earlier. 

An  advantage  of  edge  length  as  a  measure  is  that  the  penalty  is  rotationally 
invariant,  i.e.  remains  constant  throughout  all  rotations  of  the  scene  within  the  region. 
This  could  not  be  obtained  using  discrepant  pairs  as  a  measure  although  it  has  been 
shown  by  our  colleague  Robin  Sibson  that  this  variability  can  be  minimised  using  a 
down  weighting  of  \l'^2  for  the  diagonal  adjacencies. 

3.  The  Restoration  Algorithm 

The  restoration  is  done  in  three  stages,  the  first  two  ot  which  have  already  been 
described  : 

Stage  1  :  ICM  to  convergence  on  full  size  pixel  grid. 


Stage  2  :  ICM  to  convergence  on  2x2  pixel  grid. 

Stage  3  :  Updating  process  on  the  line  segments  representing  the  edges. 

Stage  3  requires  that  we  now  regard  the  reconstruction  as  a  series  of  line  segments 
separating  the  two  colours.  An  initial  representation  is  obtained  in  a  straightforward  way 
from  the  end  product  of  Stage  2.  The  updating  process  treats  pixels  in  pairs,  selecting 
the  best  place  for  two  edges  to  meet,  given  the  current  restoration  of  neighbouring 
pixels. 

As  an  example,  consider  the  configuration  at  pixels  <  and  j  shown  in  Figure  1.  The 
distances  a  and  b  are  determined  by  the  current  colouring  of  neighbouring  pixels  and 
treated  as  constant  for  the  moment  The  distance  W  is  chosen  to  minimise  the 
contribution  from  pixels  i  and  j  to  the  total  penalty  (2.1),  i.e. 


PenJWy v-  — j  £  (yk  -  Ptw)2  +  P'(«jw  +  ejw)>  (3.1) 

- 2<r  t=i,j 

where  ekW  is  the  length  of  edge  in  pixel  k  when  the  join  is  at  W  and  is  the 
proportion  of  black  in  pixel  k  when  the  join  is  at  W. 

For  the  case  shown  in  Figure  1,  this  penalty  is 

g(W)  =  -±-^yi-a-\{W-a))2Hyrb-\{V/-b))2} 


+  Q'{  v  1  +(\V-a)2+y  l+(  W-fr)2} 


Fig.  1.  Updating  the  position  of  edges  in  pixels  i  and  j. 


This  can  not  be  minimised  directly  but  the  form  of 


dg(W) 

dW 


——(IW+a-ly^b-lyj) 
4c r 


(vy-a)  .  w-b) 

Vl  +(W-a)2  'H+iW-b)2 


suggests  an  iterative  approach.  Given  an  approximate  solution  we  solve 


-+-(2Ws+a-2yi+b-2yi)  +  P' 


(W^-a)  |  (WV-fc) 

Vi +(vy,_,-a)2 


=  0 


to  obtain 


4o2p' 

a 
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+  (2y,-a+2y;-b) 

Vi+<*y,- 

i-a)2  yli+W^-bf 

2+4<t2P' 

1 

1 

Vl+(*Vi  -a)2  Vl+fHO-i-*)2 

Starting  from  any  sensible  initial  value,  W0,  accuracy  to  3  decimal  places  was 
achieved  after  at  most  four  iterations.  In  practice  we  take  W0  to  be  the  value  of  W  prior 
to  this  update. 

Different  forms  of  (3. 1)  are  possible  depending  on  which  neighbours  of  pixels  i  and 
j  contain  both  colours.  There  ate  only  four  distinct  cases  that  may  arise  and  these  are 
shown  in  Figure  2. 


Fig.  2.  Possible  configurations  of  edges  in  two  neighbouring  pixels. 


We  have  shown  the  method  of  solution  for  case  (i)  and  cases  (ii)  -  (iv)  are  solved 
in  a  similar  way.  All  other  cases  can  be  reduced  to  one  of  the  above  by  means  of 
exchanging  and/or  inverting  the  pixels  and  their  colours. 

The  most  natural  order  of  updating  the  edge  pixels  would  seem  to  be  to  follow  an 
edge  around,  updating  each  join  in  turn,  completing  circuits  of  the  edge  until 
convergence.  An  alternative  method  is  to  update  every  k,h  join  around  the  circuit, 
therefore  completing  k  laps  before  each  pixel  has  been  updated  once.  Initial  results 
suggest  that  this  provides  additional  stability  in  the  updating  process;  we  have  found  the 
value  k  =  3  to  give  particularly  good  results. 

4.  An  example 

We  illustrate^  the  methods  we  have  described  with  an  artificial  example.  Figure  3a 
shows  a  true  image  and  the  superimposed  pixel  grid.  The  record  from  which  a  restored 
image  was  was  constructed  obtained  by  generating  a  Gaussian  random  variable  for  each 
pixel  with  mean  equal  to  the  proportion  of  the  pixel  coloured  black  in  the  true  image  and 
variance  0. 12.  Figure  3b  is  the  reconstruction  after  stage  1,  in  which  the  ICM  method 
with  (3=1  has  been  used,  treating  each  pixel  as  either  completely  black  or  completely 
white.  Note  that  this  is  a  rather  poor  approximation  to  the  true  image  but  it  is  the  best 
that  can  be  done  without  dividing  pixels.  Subdividing  each  pixel  into  four  in  stage  2 
produces  the  reconstruction  in  Figure  3c:  the  amounts  of  black  in  each  full  pixel  arc  now 
much  closer  to  the  corresponding  records  and  the  divisions  of  split  pixels  match  up  well 
with  the  true  image.  Proceeding  to  stage  3,  using  (3'=2,  gives  the  final  reconstruction 
shown  in  Figure  3d,  despite  the  coarseness  of  the  original  pixel  grid  and  the  addition  of 
noise  to  the  record,  this  reconstruction  is  barely  distinguishable  from  the  true  image. 
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Fig  3a  True  image 


Fig  3b  Reconstruction  after  stage  1 
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Fig  3c  Reconstruction  after  stage  2 


Fig  3d  Final  reconstruction 


5.  Further  extensions 


(a)  Consider  a  pixel  which  has  true  colouring  as  shown  in  Figure  4.  Clearly  the 
straight  line  approximation  to  this  edge  will  be  poor  and  could  have  an  adverse  effect  on 
the  reconstruction  of  neighbouring  pixels  and  pixels  further  along  the  edge.  This  may  be 
overcome  using  a  more  intricate  restoration  method,  e.g.  allowing  two  straight  lines 
meeting  at  some  point  within  a  pixel. 
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Fig.  4.  A  pixel  containing  a  boundary  that  can  not  be  approximated  well 
by  a  single  straight  line. 


(b)  The  method  presented  in  this  paper  can  be  extended  to  scenes  containing  more 
than  two  different  colours.  Where  any  two  regions  meet  we  can  adjust  the  algorithm  to 
provide  a  continuous  line  join.  More  computation  is  required  to  find  the  best  colouring 
for  a  pixel  in  whicn  three  or  more  regions  meet 
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